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Three-Dimensional  Seismic  Ray  Tracing  in  a  Laterally 
Heterogeneous  Spherical  Earth1 

Kiju  n  H.  Jacob 

Lnmont-Dolierly  (leitlnyical  Observatory  of  Columbia  University 
Palisades,  Sew  York  W!Mi/t 

Recent  seisniological  studios  suggest  lutrriil  inhomogencities  in  P  and  S  velocities  of  the 
mantle  that  are  associated  with  slabs  of  mobile  lithosphere  descending  into  tlu  mantle  be¬ 
neath  island  ares.  In  special  eases,  travel  times  of  P  traversing  such  zone's  can  differ  bv  as 
much  as  5  see  and  of  S  by  up  to  10  see  from  standard  travel  times.  In  addition,  such  zones 
are.  characterized  by  relatively  low  attenuation  of  iS’-wave  energy  compared  with  high  attenua¬ 
tion  in  a  broad  zone  on  the  landw’ard  side  of  the  active  volcanoes.  To  explain  the  observed 
anomalous  travel  times  and  attenuation  phenomena,  it  is  necessary  to  trace  the  path  of  body 
waves  through  laterally  heterogeneous  earth  models.  The  technique  of  ray  tracing  developed 
here  uses  Fermat’s  principle  to  obtain  the  differential  equation  of  a  ray  in  spherical  coordi¬ 
nates.  The  position,  direction,  and  travel  time  of  the  seismic  wave  front  at  any  point  along 
the  curved  ray  path  are  obtained  by  numerical  integration  of  the  differential  equation  for  an 
assumed  three-dimensional,  continuous  velocity  distribution.  The  problem  of  representing  a 
realistic  three-dimensional  velocity  structure  in  the  earth  is  solved  in  a  way  that  is  especially 
suitable  for  use  on  computers.  Some  examples  for  rays  traversing  an  island-arc  structure  are 
presented.  The  implications  of  this  method  of  tracing  rays  in  a  laterally  heterogeneous  earth 
are  discussed  with  respect,  to  seismic  travel-time  studies,  interpretation  of  residuals  in  terms 
of  tectonic  heterogeneities,  source  bias,  and  the  precise1  location  of  earthquakes  and  nuclear 
explosions;  tlT/tll  measurements  from  large  seismic  arrays  and  their  inversion  to  obtain  de¬ 
tails  of  the  velocity  structure  in  the  upper  mantle  are  also  discussed. 


Since  the  beginning  of  instrumental  seismol¬ 
ogy  and  throughout  the  first  half  of  this  cen¬ 
tury,  seismologists  have  usually  treated  the 
seismic  velocity  structure  of  the  earth’s  interior 
as  spherically  symmetric.  It  was  not  until  1P33 
that  fSutenborg  suggested  that  errors  due  to  the 
earth’s  elliptirity  might  he  significant  for  telc- 
seismie  travel  times,  mid  hence  cllipt icily  should 
ihen  1h>  taken  into  account  | Gutenberg  and 
Richter,  1!KM).  Jeffreys  subsequently  derived  a 
method  to  correct  ’ravel  times  for  clliptirity 
for  any  given  epicenter-station  configuration 
\JefJniis,  I  MSI.  lie- all's  cllipt  tcity,  no  other 
lateral  variations  in  the  seismic  velocities  of  the 
up)MT  mantle  were  emphasized  until  10  or  JO 
years  ago,  though  regional  variations  iu  the 
thickness  and  structure  of  t lie*  continental  and 
oceanic  eru-t  gradually  Inca  me  evident 

Numerous  geophysical  stndie-  iu  the  la-t  two 
decade-  demon-irate  a  eo!l«i*tenl  global  put  - 

1  l.auiol|l-l)olii  rty  <  «<  olo||0  ;d  (Nm-rvaiorv  Col<- 
inbmioii  1580 

I  ••(<%  it  iflii  c  1  «:*«  In  iht  .\ni«  iu. <ft  <  t'lipli)  .«!  1  mod 


tern  of  strong  lateral  variations  of  seismic  ve¬ 
locities  and  other  physical  projierties  not  only 
in  the  earth’s  crust  but  also  of  the  upper  700 
km  of  the  mantle.  These  heterogeneities  are 
closely  related  to  major  tectonic  features  such 
as  island  ares,  mid-oceanic  ridges,  rift  and  frac¬ 
ture  zones,  and  orogenic  Is'lts.  It  is  surprising 
that,  despite  the  accumulating  evidence,  no 
suitable  method  has  1h>cii  available  until  now 
for  calculating  ray  paths  ami  travel  times  of 
seismic  ImhIv  waves  through  a  laterally  hetero¬ 
geneous  earth. 

I*  is  the  purpose  of  this  -indy  to  fill  that 
gap  bv  presenting  a  basic,  computer-oriented 
method  of  tracing  seisin'.-  rays,  either  for  P 
er  .S  waves,  propagating  through  a  laterally 
heterogeneous  spherical  earth.  Tin-  method  i- 
csjiecudly  designed  to  consider  a  class  <<1  lateral 

-ci-mic  lieterogeucil  h  s  r  lated  . . a-tloor 

•pleading  ami  global  tec  nine-,  eg  ,  the  deep 
-truelure  In 'Heath  lllld-oeealuc  ridge-  ami  l-- 
l.i i o  1  ar.  -.  i be  -oiiree-  and  -ink-,  re-pi'iivelv,  m 

a  -\-lem  ol  dnt’ll  g  ll'lio-phene  plalcv 

Th«  plan  lor  tin-  pa|sr  i-  (It  to  de-erils' 
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t ho  geometric  representation  of  lateral  hetero¬ 
geneities;  (2)  to  derive  tlio  differential  equa¬ 
tions  for  a  ray  from  Fermat’s  principle ;  (3)  to 
present  a  numerical  scheme  for  the  integration; 
(4)  to  disenss  some  of  the  ray  characteristics  in 
a  laterally  heterogeneous  earth;  (5)  to  give  an 
exa  mple  for  applying  the  ray-tracing  method  to 
travel-time  studies;  and  (6)  to  disenss  some 
of  the  seisinologirai  implications  of  the  new  ray¬ 
tracing  method. 

Geometrical  Representation  of  Lateral 
Heterooeneities 

Any  point  in  a  spherical  earth  can  lie  described 
by  means  of  a  spherical  coordinate  system, 
r,  6,  \;  r  is  the  distance  (in  kilometers)  from 
the  earth’s  center,  6  is  the  polar  distance  (in 
degrees)  from  the  north  pole,  and  A  is  the 
geographical  longitude  (in  degrees),  with  posi¬ 
tive  A  east  of  Greenwich  and  negative  A  west 
of  Greenwich.  The  limits  for  the  coordinates 
r,  6,  A  are 

0  km  <  r  <  R 

0°  <  6  <  180°  (1) 

-180°  <  A  <  +180° 

R  is  the  radius  of  the  earth. 

Major  lateral  heterogeneities  of  the  earth  are 
associated  with  the  tectonic  structure  of  island 
arcs  and  mid-oceanic  ridges.  The  seismicity  and 
tectonic  features  of  island  arcs  are  described 
in  detail  by  Katsumata  [I960],  Sykes  [1966], 
Oliver  and  Isacks  [1967],  Vtsu  [1967],  and 
Mitronovax  el  at.  [1969]  among  others.  Tlx* 
inclined  zone  of  seismic  foci  approximately 
coincides  with  the  descending  lithospheric  plate 
that  is  characterized  bv  low  attenuation  (high 
Q)  and  high  seismic  velocities.  Such  slabs  of 
lithosphere  extend  from  the  trench  down  into 
the  mantle  Itcncath  the  licit  of  active  volcanoes 
to  a  maximum  depth  of  alsmt  700  km.  In  con¬ 
trast  to  the  descending  plate,  low  seismic  mantle 
velocities  and  high  a  Henna  lion  (low  Q)  are 
ob*erved  in  a  broad  zone  on  the  landward  side 
el'  the  \oleamc  Iieli. 

To  map  i lie  subsurface  Imnndaric*  of  lie 
different  i  eel  nine  uniis,  n  is  convenient  io 
delenmne  contour  lilies  :ii  varioii'  depth* 
h\H.  Al  _  cofl'tutll  Oil  I  lie  interface*  -eparallllg 
nljaeenl  belerogclleoll*  Iiodle*  \*  ]Hiilllei|  oil* 


before,  the  seismic  zone  beneath  an  island  arc 
approximately  defines  the  dipping  lithospheric 
plate.  Sykes  [I960],  Sykes  ct  at.  [1969],  and 
Katsumata  ami  Sykes  [  1969]  projected  precisely 
relocated  hypocentcrs  on  vertical  sections  per¬ 
pendicular  to  that  seismic  zone  in  several  island 
arcs  in  the  Pacific.  Similar  studies  were  carried 
out  by,  among  others,  Hamilton  and  dale 
[1968]  for  New  Zealand  and  Katsumata  [1960] 
and  I'tsu  [1967]  for  Japan.  From  such  sections 
or  from  hypocenter  maps  it  is  possible  to 
contour  the  surface  of  the  dipping  seismic  zone 
in  each  active  arc.  Figure  1  shows  as  an  example 
a  contour  map  of  the  seismie  zone  in  the  Fiji— 
Tonga-Kermadee  arc.  The  contours  were  de¬ 
rived  from  seismic  events  reported  by  the  U.S. 
Coast  and  Geodetic  Survey  (COS)  for  the 
period  1961  to  1968.  To  store  information  on 
the  depth  contours  for  use  in  a  computer,  longi¬ 
tudes  Ai ,  are  read  for  the  intersections  of  each 
contour  h,  with  equally  spaced  colatitndes  6,. 
This  sampling  procedure  is  suitable  only  for  pre¬ 
dominantly  N-S  striking  tectonic  features.  For 
an  K-W  striking  feature,  a  corresponding 
scheme  has  to  be  used;  polar  distances  0 ,,  are 
determined  for  the  intersections  of  contours  h, 
with  equally  spaced  longitudes  A,. 

Table  1  contains  the  longitudes  Af,  sampled 
along  the  contours  It,  for  the  seismie  zone  in  the 
Fiji-Tonga-Kermadec  area  shown  in  Figure  1. 
This  table  demonstrates  schematically  how  the 
information  of  an  arcuate  structure  is  stored  in 
the  computer.  So  far,  however,  only  an  arcuate 
dipping  surface  has  lieen  represented.  To  pro¬ 
vide  the  dipping  plate  with  some  volume,  one 
has  to  assign  horizontal  thicknesses  a  and  b, 
resjiectivcly,  on  each  side  of  tin*  seismic  surface, 
which  will  Is*  labeled  .8  from  now  on.  This 
configuration  is  schematically  deni  oust  rated  in 
Figure  2  with  a  section  of  an  arc.  The  hori¬ 
zontal  distances  a  and  b  are  measured  |x*r]icii- 
dicular  io  the  strike  of  each  depth  contour  h, 
and  for  each  colaiitnde  interval  1 0,,  0,.,\.  The 
thickness  a  is  on  the  concave  side  of  the  seismic 
surface  (toward  the  volcanoes  I,  and  b  is  on  ihe 
convex  side  (toward  lhe  ocean)  The  Iwo 
defined  liound:irie»  of  the  dipping  plate  are 
called  .1  and  H,  re*|iecii\oly ;  ihe  individual 
facci.  conipri'iiig  .1  ami  H  are  parallel  lo  l lie 
corn-* | winding  facei*  ot  8  in  each  depth  interval 
|  h..  h  |  and  in  each  jwilar  dl»iaiicc  interval 

!«.«.,  I 
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Kite.  1.  Representation  of  (lie  dipping  seismic  zoni>  in  the  Tonga-Kerinadec  island  arc  by 
means  of  depth  contours  h  =  constant.  Tlic  numbers  on  the  contours  arc  depths  in  kilometers 
below  the  earth’s  surface.  Open  circles  show  points  of  intersection  with  latitude  <p  —  constant ; 
their  longitudes  \  are  compiled  in  Table  1. 

In  the  unlikely  case  that  the  plate  thickness  imposed  on  a  standard  velocity  distribution  for 

or  distance  bet  ween  S  and  .4  or  II  is  of  the  same  the  mantle  |c.g.,  Jeffreys  ami  Hullen,  11)40, 

order  of  magnitude  ax  the  radius  of  curvature  ISMS;  Ilerrin  et  a!.,  ISMiSj,  which  is  a  continuous 

of  the  arc,  then  A  or  II,  respectively,  is  no  function  of  depth  only. 

longer  a  smooth  boundary  but  is  an  irregular  The  existence  of  a  region  of  low  seismic 
surface  composed  of  facets  with  small  offsets  at  velocities  beneath  and  adjacent  to  the  volcanic 
each  sampling  value  of  the  polar  distance  0,.  licit  was  discussed  by  Mitrunuvas  |T!)fl!)|,  I'tsu 

These  offsets  can  lx*  kept  small,  however,  if  a  |11M>7|,  M ulnar  am!  Oliver  |11M>1I|,  and  Kana- 

xnflicicntly  small  sampling  interval  <16  is  used.  won  |  ID70|,  among  others.  Karig  |1!)7()|  found 

A  certain  degree  of  irregularity  in  the  interfaces  indications  of  young  extcnxional  tectonic  fea- 

.4  and  II  can  lie  tolerated  Ijoeause  the  offsets  turcs  to  the  west  of  the  volcanic  belt  in  the 

are  smoothed  in  the  specific  numerical  approach  Tonga-Kcrmadcc  and  Philippine  arcs,  Thus 

used  for  describing  individual  seismic  rays.  For  the  zone  of  low  velocity  behind  island  arcs 

the  dipping  plate  limited  by  the  surfaces  ,4  might  l>c  related  to  some  kind  of  local  ocean- 

anil  II,  a  seismic  velocity  is  assigned  that  is  door  spreading.  To  account  for  such  a  region  of 

higher  by  a  certain  |H-reentage  &rtl  than  the  lower  velocities,  another  interface  < '  is  intro- 

veloeity  ai  the  corresponding  depth  in  the  duecd  to  form  a  wedge-shaped  body  of  depth  </, 

mantle  outside  the  plate.  Hence  the  boundaries  and  of  width  r„  at  the  earth’s  surface,  as  shown 

.4  and  H  represent  discontinuities  in  seismic  in  Figure  ‘J.  The  facets  of  the  interface  (‘  al>o 

velocities.  Tlii-  velocity  differential  is  sii|»er-  parallel  the  seismic  zone  S  in  the  interval' 
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TABLE  1.  Numerical  Represent  at  ion  of  the  Seismic  Zone  of  the  Tonga-Kermadec  Island  Arc 


Longitudes,  deg 


South 

Latitude, 

deg 

h  - 

700  km 

h  = 

600  km 

h  - 
500  km 

h  - 

400  km 

h  - 

300  km 

h  = 

200  km 

h  - 

100  km 

h  = 

0  km 

14 

•  •  • 

•  •  • 

•  •  • 

,  ,  , 

.  •  • 

.  .  . 

-174.8 

-173.2 

15 

.  .  . 

.  .  . 

.  .  . 

.  .  . 

— 177.3 

-175.0 

-174.0 

-172.5 

16 

.  .  . 

.  .  . 

.  .  . 

-179.0 

-176.9 

-174.8 

-173.6 

-172.4 

17 

-179.5 

-170.0 

-  178.7 

-179.0 

-177.2 

-174.5 

-173.6 

-172.5 

IS 

-178.0 

-178.6 

-178.1 

-178.0 

-176.7 

-174.9 

— 174.0 

-172.8 

10 

-178.7 

-178.3 

-177.5 

-177.5 

-176.5 

-175.1 

-174.5 

-173.1 

20 

-179.0 

-178.5 

-177.5 

-177.5 

-176.5 

-175.7 

-175.0 

-173  5 

21 

-179.5 

-179.1 

-178.3 

-177.5 

-176.5 

-176.1 

-175.4 

-173.0 

22 

+  179.7 

-170.8 

-179.4 

-178.2 

— 177.0 

-176.6 

-175.8 

-174.3 

2d 

+  179  1 

+  179.7 

-179.7 

-178.8 

— 177.7 

-177.1 

-176.2 

-174.6 

24 

+  178.6 

+  179.3 

-179.9 

-179.3 

-178.2 

-177.4 

-176.3 

-175.0 

25 

+ 178.3 

+  179.1 

+  179.9 

-179.4 

-178.5 

-177.5 

-176.4 

-175.2 

26 

+  178.1 

+  178.7 

+  179.5 

—  179.5 

-178.8 

-177.7 

-177.5 

-175.3 

27 

+  178.4 

+  179.2 

-179.5 

-178.9 

-177.0 

-177.5 

-175.6 

28 

+  179.1 

-179.5 

-178.9 

-178  1 

-177.6 

-175.8 

20 

-179.5 

-178.9 

-178.4 

—  177.7 

-176.3 

:io 

-179.6 

-179.2 

-178.7 

-178.0 

-176.8 

dl 

-179.9 

-179.5 

-179.0 

-178.3 

-177.2 

d2 

+  179.8 

-179.7 

-179.3 

-178.8 

—  177.6 

dd 

+  179.5 

+  180.0 

-179.6 

-179.1 

-178.1 

34 

+  178.8 

+  179.6 

+  180.0 

-179.1 

-178.5 

d.r> 

+  178.4 

+  178.7 

+  179.3 

+  179.7 

-179.0 

d6 

+  177.6 

+  178.2 

+  178.6 

-179.8 

d7 

+  176  6 

+  177.2 

+  177.9 

+  179.7 

d8 

+  175.6 

+  176.2 

+  177.0 

+  179.1 

d'J 

.  .  . 

+  175.0 

+176.1 

+  178.3 

40 

.  .  . 

.  .  . 

+  175.2 

+  177.4 

Longitudes  are  sampled  at  the  intersections  of  depth  contours  h  «  constant  with  parallels  at  latitudes 
<t>  =  constant. 


\0,,  The  horizontal  distance  between  S 

and  ('  at  any  depth  h  <  d  is  defined  as 

c(h)  =  c0  —  (cn  —  a)h/d  (2) 

and  is  measured  perpendicular  to  the  local 
strike  of  S.  According  to  (2),  C  approaches  A 
at  depth  d,  indicated  by  the  point  D  in  Fig¬ 
ure  2.  The  wedge-shaped  volume  is  bounded  by 
A,  and  the  surface  of  the  earth.  The  velocity 
inside  this  wedge  can  be  chosen  to  be  a  certain 
percentage  Stv  lower  than  the  normal  mantle 
velocities  for  any  depth  h  <  d.  To  describe  the 
gross  seismotectonic  features  of  an  island-arc 
system  as  presently  conceived,  it  is  thus  neces¬ 
sary  to  define  only  six  parameters  a,  b,  c„,  d, 
St',,.  and  « tv.  In  addition  the  coordinates  of 
the  depth  contours  of  the  seismic  zone  S  must 
1h>  specified. 


Fermat’s  Principle  and  Euler’s 
Differential  Equations  of  the  Ray 

The  main  purpose  of  this  study  is  to  provide 
a  method  of  tracing  seismic  rays  through  an 
arbitrary  velocity  structure.  To  some  extent  we 
follow  the  approach  of  Sattlegger  [1964],  who 
treated  a  similar  problem  in  a  Cartesian  space 
for  seismic  reflection  surveys.  Instead  of  using 
the  velocity  function  v(r,  6,  A)  to  describe 
the  velocity  structure  of  the  earth,  we  use 
here  the  more  suitable  quantity  slowness 
u(r,  6,  A),  the  reciprocal  of  the  velocity.  We 
assume  that  the  scalar  u  is  a  continuous  func¬ 
tion  throughout  the  earth.  Hence,  any  first-order 
discontinuity  in  slowness  u  is  replaced  by  a 
zone  of  gradual  change  in  u.  To  allow  applica¬ 
tion  of  ordinary  rav  theory,  thn  magnitude  of 
the  gradient  of  u  must  be  kept  reasonably  small. 
Thus,  for  a  given  contrast  in  slowness  u,  the 


KKIMMIC  HAY  THACINO 


6679 


zone  of  gradual  change  in  u  lias  to  ho  wide 
enough  to  keep  the*  gradient  of  u  small. 

According  to  Fermat's  principle  a  seismic 
ray  will  elmose  the  specific  path  S  that  takes 
the  least  travel  time  T  to  propagate  from  a 
point  A  to  a  point  B : 

CH 

T=  ud.I.Y  3  ds  =  Extremum  (3) 
J  A 

where  /,  is  the  tangent  vector  at  any  point 
along  the  ray  with  its  components 

/,  =  (r',  rO',  r  sin  0\')  (4) 

Tfie  tph uiu ir«  f'  iV'  '  ’m>  'he  drfiVuliui  ut 
the  coordinates  along  the  ray  increment  ds  and 
represent  the  components  of  direction  of  the 


ray.  Multiplied  by  their  respective  metric 
factors  I,  r,  and  r  sin  0  for  a  spherical  coordi¬ 
nate*  system,  they  form  the  components  of  a 
unit  tangent  vector  where 

(M,)1/2  =  1  (5) 

and  the  summation  convention  applies. 

The  problem  stated  in  equation  3  is  common 
in  the  calculus  of  variation.  It  can  be  solved 
by  finding  solutions  to  the  corresponding  set  of 
Euler’s  differential  equations.  They  are  obtained 
by  defining  the  integrand  «(/,/, ),/J  of  equation  3 
as  the  Eider  function  F(g,,  g,'),  where  g,  stands 
fur  iUm  gi'uonlinil  r uordimTrM  itwt  f kt  rh«r 
derivatives;  the  differential  equations  take  the 
following  form: 


big.  2.  Section  of  an  island  arc  structure  showing  schematically  the  geometry  of  the  model 
used  for  representing  lateral  heterogeneities  related  to  the  deep  seismic  zone.  S  =  seismic 
zone ;  A  —  upper  boundary  of  the  descending  plate  of  lithosphere;  B  ~  lower  interface  of  the 
descending  plate;  (’  =  boundary  on  the  continental  side  of  the  wedge-shaped  low-velocity 
body  near  the  volcanic  belt;  I)  ~  decj  cst  point  of  the  low-velocity  wedge  at  depth  d.  Open 
circles  on  .S  represent  sampled  points  shown  also  in  Figure  1  and  Table  1.  For  details  of  meas¬ 
uring  tin'  horizontal  distances  a,  h,  and  e.  see  text.  The  (|iiestion  mark  in  the  upper  part  of 
B  indicates  that  the  model  does  not  represent  here  the  actual  boundary  of  the  lithospheric 
plate  very  well. 
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K, 


d 

(In 


From  (-1),  (ft),  ami  (6)  wo  obtain 


(0) 


T  («'•')  -  ?-  -  >"‘(0"  +  A'2  «i"2  9)  =  0 
da  dr 

y  («r20')  —  ~  —  JurV2  sin  20  =  0  (7) 

da  d  u 


y  («r2X'  sin2  0)  —  ^  =  0 
as  oX 

Tlio  roiiiaiuiiiK  differentiation  of  flic  first  Icnn 
ami  rearranging  of  all  resulting  terms  yields  the 
following  expressions  for  I  Ik*  second  derivatives 
r”,  ft",  X": 

r"  -  0/«)(fj!  -  »v) 

+  r(0'2  +  X'2  sin2  6) 


0"  =  -i  ~  +  *X'2  sin  20 
ur  06 


-Kr  +  ’r) 


(H) 


X"  = 


1  du 


ur2  sin2  0  OX 


,/u'  r'  O'  \ 

-  X'  -  +  2  L  +  2 

\w  r  tun  0/ 


where 


,  du  du,  I  du  .  1  du  , 

u  ~  a  —  A  r  ^ —  y,  0  H - :  o  y  A 

os  or  r  dO  r  sin  0  dX 


(9) 


The  (‘(plat ions  s  are  the  differential  equations  of 
the  seismic  ray  in  terms  of  components  of  ray 
curvature  r" ,  0",  X"  anywhere  in  a  spherical 
earth  that  has  a  continuous  slowness  distribu¬ 
tion  nlr,  0,  X).  The  ray  curvature  depends, 
first,  on  the  components  of  the  gradient  of 
slowness  (Am  dr),  (du  00),  (Ok  OX)  and,  second, 
on  the  direction  of  the  ray  with  respect  to  the 
direction  of  the  gradient  of  slowness.  The  second 
relation  is  concealed  in  equation  0  for  it'  = 
dn  ds.  It  represents  the  projection  (Scalar  pro- 
duel  )  of  the  gradient  of  slowness  on  the  tangent 
vector  of  t lie  ray. 


K.X  CAN  SION  OK  H.w  (  00111)1  NATHS  INlo  A 
Tavuiii  Skuiks 

Any  attempt  to  find  the  actual  ray  path  l>\ 
analytical  integration  of  (N )  for  a  realistic 
slowness  function  u(r,  0,  X)  is  likely  to  fail. 
Instead  of  an  analytical  integration  of  (N),  we 
try  stepwise  numerical  integration  by  a  finite- 
series  scheme,  a  method  similar  to  that  used  by 
Sat  tinnier  |  i‘ili-1,  I  Olid  |. 

The  behavior  of  a  seismic  ray  s  in  the  vicinity 
of  a  point  0„,  X„)  on  the  ray  can  be 

described  by  expanding  the  ray  coordinates 
r,  0,  X  and  travel  time  /  in  a  Taylor  series: 

r  =  r„  +  r,/  da  -f  W  da*  +  •  •  • 

0  =  0.i  H-  0>/  ds  +  j0„"  ds‘  +  •  •  • 

X  =  X„  +  X,/  da  +  JX,,"  da*  +  •  •  • 

/  =  /„  +  m„  da  +  }tu»  c/s2  + 

All  values  forming  the  zero-,  first-,  and  second- 
order  terms  on  the  right-hand  side  of  (10)  are 
either  known  or  can  be  calculated  from  (N)  and 
(0)  if  th(>  slowness  u(r,  0 ,  X)  is  given  ;  /•„,  0„,  A,„ 
and  /„  represent  the  initial  position  and  the 
origin  time  of  the  ray  at  the  point  /V;  r.f,  0,,',  X,,' 
represent  the  initial  direction  of  the  ray  that 
can  be  prescribed.  Finally,  0„",  X„"  are  the 
elements  of  ray  curvature  and  can  be  calcu¬ 
lated  from  (S)  and  (it).  If  the  spatial  gradient 
of  slowness  //  varies  slowly,  the  higher  order 
terms  do  not  contribute  significantly  to  the  sum 
in  (10).  The  error  in  neglecting  these  terms 
call  be  kept  arbitrarily  small  by  proceeding 
only  in  small  increments  ds.  After  the  first  step 
da,  the  ray  lias  propagated  from  /•*„(/•,,,  0„,  X„,  t„) 
to  l‘(r,  0,  X,  t).  The  new  direction  of  the  ray 
point  /'(r,  0,  X,  t)  can  be  obtained  by  differ¬ 
entiating  (10); 

r'  =  r,/  +  r„"  da  +  •  •  • 

O'  =  0,/  +  ft,"  da  +  •••  (11) 

X'  =  X,/  +  X„"  da  +  •  •  . 

Only  the  tir>t  two  terms  in  (II)  must  be  con¬ 
sidered  to  keep  the  accuracy  of  the  computa¬ 
tion  consistent  with  that  in  equation  10.  Terms 
of  higher  order  can  be  neglected  again  for  the 
proper  choice  of  da. 

I  'sing  r,  I),  X  from  (10),  r',  O',  X’  from  (II), 
and  /",  0”,  X"  from  <M,  we  may  trace  the  ray 
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from  point  to  point  in  steps  ds  until  the  ray 
emerges  somewhere  at  the  earth's  free  surface. 

Numerical  Scheme 

To  let  the  ray  propagate  point  to  point,  say 
from  Pk  to  PA.,,  one  has  to  determine  the  compo¬ 
nents  of  eurvature  of  the  ray  r",  9",  X",  or, 
aeeording  to  equation  S,  to  find  tlte  gradient 
of  slowness  ( du  dr,  du'dO,  du  <DA)  in  the  vicin¬ 
ity  of  Pi.  It  is  obvious  that  the  gradient  of  u  at 
an  intermediate  point  Pk.„t  :s  more  representa¬ 
tive  ot  the  eurvature  for  the  ray  path  between 
Pk  and  Pi. i  than  the  gradient  at  the  end  point 
Pi.  Thus  the  aeeuraey  of  the  numerieal  integra¬ 
tion  could  he  improved.  Unfortunately,  it  is 
not  possible  to  determine  the  loeation  of  the 
point  Pi.,  ■  and  consequently  the  gradient  of 
this  point,  because  the  eurvature  for  the  ra\ 
between  I\  and  Pi.,  is  not  known  in  advance. 
To  find  an  approximate  location  of  Pt.,,.  how¬ 
ever,  we  may  use  the  curvature  of  the  previous 
ray  element  determined  at  Pk.,,j  and  extrapo’ate 
the  ray  path  beyond  I\  by  half  an  increment 
(is/ 2.  Accordingly,  the  following  numerical 
scheme  is  used  for  the  ray-tracing  proeedure: 

1.  Find  the  auxiliary  point  Pi.,  a  by  extrap¬ 
olating  the  previous  ray  element  beyond  Pk  by 
an  increment  c/s.  2: 

r*  +  )/j  =  rk  +  \r/  ds  +  \rk.,,/'  ds 2 

6k* i.2  =  6k  +  \6/  ds  +  \6k.,  /'  ds 2  (12) 

1/2  =  'Xi  +  ^X/  ds  +  jXt-i  i"  ds 1 

2.  Determine  the  slowness  of  u  and  its  gradi¬ 
ent  {du  dr,  du  d6,  du  dX)  at  the  auxiliary  point 
Pi., ,  and  use  these  data  in  equations  S  and  0 
to  calculate  the  curvature  components. 

3.  Let  the  ray  propagate  from  Pk  to  Pk., 
with  a  step  ds  according  to 

G* i  =  rk  +  r/  ds  +  §r*.,  "  ds 2 

6k.\  =  6k  +  6/  ds  +  \6k  + 1  i"  ds 2 

X* . ,  =  Xt  -j-  X/  ds  +  jX,.i  /'  ds 2 

/*.  i  =  h  +  a*. i/2  ds 

4.  The  direction  of  the  ray  at  the  new  point 
Pi.,  i-  given  by 

’•*./  =  rk  +  r*.I  2"  ds 

6k./  =  ek'  +  ek.l  2"  ds  (14) 

X*./  =  X/  +  X*.,  ds 


5.  Repeat  steps  1  through  4  with  all  sub¬ 
scripts  increased  by  one  integer  unit,  and  con¬ 
tinue  until  the  ray  reaches  the  earth’s  free 
surface. 

The  slowness  and  its  gradient  components  are 
determined  for  step  2  of  that  scheme  by  one 
of  the  following  two  methods,  depending  on 
whether  the  ray  penetrates  a  zone  with  or 
without  lateral  heterogeneities  present. 

(a)  In  the  presence  of  lateral  heterogenei¬ 
ties:  Assume  a  cube  centered  around  the 
auxiliary  point  Pk.,  .  with  an  edge  length  2ie. 
Its  eight  corner  points  are  labeled  Q,m„  where 
the  -uhseript  /,  m,  or  n  may  take  either  the 
value  1  or  2,  depending  on  whether  its  corre¬ 
sponding  coordinate  r,  6,  and  A  of  the  point 
Q,„„  is  smaller  or  larger  than  that  of  the  auxili¬ 
ary  point  Pi.,  i. 

Thus  these  eight  points  Q,,„„  are  located  at 

Qm  =  (r  —  u\  6  —  w/r,  X  —  u>/r  sin  6) 

Q211  =  (r  +  tc,  6  —  w/r,  X  —  w/r  sin  0) 

Qm  —  (r  +  u\  6  +  w/r,  X  +  w/r  sin  6) 

(15) 

The  slowness  at  each  point  Qlm„  is  called  u 
From  these  eight  values  of  slowness  we  derive 
the  arithmetic  mean  and  assign  it  to  the  slow¬ 
ness  at  point  Pi.,, 2. 

«*.l  2  =  (m hi  +  «112  +  Um  -f  U,22 

+  W21I  +  U212  +  U221  +  «222)/8  (16) 

By  calculating  the  average  of  the  slowness  of 
4  [mints  at  the  higher  value  of  one  specific 
coordinate  and  subtracting  it  from  the  average 
slowness  of  4  points  at  the  lower  value  of  that 
coordinate,  one  can  obtain  the  components  of 
the  slowness  gradient  for  this  coordinate: 

dll  6r  =  [(u222  +  Uj21  +  «2I2  +  «2ll) 

-  (Will  +  W,21  +  U|,2  +  U|22)]/8u1 
du  d6  =  [ ( w  1 2 1  +  w,22  +  “221  +  N222) 

—  (um  -F  w ,  1 2  ■+■  w 2 1 1  +  U212)]  (8 it'  r) 
du  d\  =  1(m,,2  +  u  1 22  +  W212  +  “222) 

-  («ni  +  w , 2 1  +  Wan  +  M22i)]/(8uv  r  siii  6) 
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II  is  now  easy  to  understand  liow  tin  slmi]' 
diseontinuities  in  slowness  u  associated  with 
the  teetonie  features  deseribed  in  a  previous 
section  are  smoothed  out  to  gradient  zones  with 
a  finite  thickness.  That  thickness  varies  between 
2ie  and  2 w('.i)va,  de|)ending  on  the  relative  local 
orientation  of  the  auxiliary  euhe  Q(„,„  that  is 
carried  along  with  the  ray  as  it  propagates. 
The  auxiliary  eube  detects  the  orientation  of 
the  interface  by  averaging  over  the  gradients 
determined  at  successive  points  /\_1/2,  Plt ,/s, 
.  .  .  ;  during  each  step  the  gradient 
depends  on  which  of  the  eight  points  Qi„,„ 
are  found  to  be  inside  a  lateral  heterogeneity 
or  outside  in  the  regular  mantle.  To  guar¬ 
antee  a  satisfactory  averaging  of  gradients  while 
passing  through  a  transition  zone,  the  ray  must 
propagate  for  many  steps  inside  that  gradient 
zone.  This  can  be  assured  if  the  step  size  ds  is 
kept  small  compared  to  the  length  2tc  of  the 
gradient  eube.  Consequently,  a  provision  has 
been  made  in  the  computer  program  to  detect 
whether  any  of  the  components  of  the  gradient 
exceeds  a  limit  hv  using  the  following  criterion: 


du 

dr 


+ 


I  — 

r  80 


+ 


_ 1  du 

r  sin  6  d\ 


<  G 


(18) 


where  G  is  of  the  order  of  0.1  x  10'4  see/knr. 

If  cond'tion  IS  bolds,  no  unusual  boundary 
is  present  and  ds  is  generally  chosen  to  lx*  about 
1  to  2  kin;  the  eube  iengtli  2 w  is  kept  *,f  the 
order  of  2  to  4  km.  If,  however,  condition 
IS  does  not  hold  because  a  strong  gradient 
zone  is  being  approached,  then  ds  will  be 
lowered  to  0.1  or  0.2  km  while  keeping  2 w  con¬ 
stant  at  2  to  4  km.  This  allows  for  at  least 
20  steps  for  ray  propagation  through  a  gradient 
zone  replacing  the  sharp  discontinuity. 

( b )  In  flic  ubsenee  of  lateral  heterogenei¬ 
ties:  In  the  ease  where  there  are  no  lateral 
heterogeneities  present,  the  calculation  of  m4.,  * 
and  of  its  gradient  can  be  reduced  considerably 
by  using  the  slowness  u..  and  «,  of  only  two 
points  just  above  and  below  the  auxiliary  point 


Qj  ~  0\  .12  +  0 1 .  i  a,  X* , ,  i) 

Q i  =  (r».i  i  m\  Ot .  ]  j,  A( .  i  2) 


(19) 


The  slowness  u  assigned  to  the  vicinity  of  the 
auxiliary  pi, tut  /', ., ..  is 


and  the  components  of  the  gradient  of  u  are 
du/dr  =  (u2  -  m,)/2u>  ^ 

du/dd  =  du/d\  =  0 

It  is  important  to  note  that  the  method  used 
here  breaks  down  in  a  narrow  cylinder  extend¬ 
ing  from  the  south  to  the  north  pole  along  the 
rotational  axis  of  the  spherical  coordinate  sys¬ 
tem  because  of  singularities  in  the  denominator 
of  some  terms  of  (S)  and  (9)  at  r  =  0  and 
0  —  0°  mid  6  =  ISO0.  Rays  which  penetrate 
into  that  cylinder  with  a  critical  radius  of 
about  10  km  must  be  omitted. 

Additional  Information  from  thk 
Ray  Trace 

The  rav-tracing  procedure  as  described  above 
provides  information  on  the  ray  r,  0,  A  at 
discrete  points  I\,  the  direction  r',  O',  A'  at 
these  points,  the  travel  time  t,  and  the  cumu¬ 
lative  paths  along  the  ray.  From  some  of  these 
data  and  the  stored  data  on  slowness  u(r,  6,  A), 
additional  information  can  be  obtained  to 
describe  some  ray  characteristics  that  are  com¬ 
monly  used  in  seismology,  e  g.,  dt/dA  along  the 
ray,  azimuth  a  (measured  clockwise  from  north 
in  a  horizontal  plane),  and  angle  of  emergence  i. 

Relations  between  U,  O',  A'  and  the  strike  a 
and  angle  01  emergence  i  (as  measured  from 
the  downward  vertical)  can  be  found  from 
simple  geometry: 

r'  —  —  cos  t 
O'  =  —  sin  t  cos  a  It 
A'  =  sin  t  sin  a  (r  sin  0) 

The  corresponding  inverse  relations  are 

a  =  tun-1  (A'  sin  0/  —  O') 

i  =  tan-'  [r( A'2  sin*  0  +  0’*)'  *; -r'} 

Note  that,  while  «  may  lie  in  any  of  the  4 
quadrant.',  i  can  lie  only  in  the  first  and  second 
quadrant . 

An  imporiam  qiianiiiy  111  cla>Mc:il  ray  I  henry 
will)  spherically  symmetric  velocity  >irnciurc 
etc  I  i>  lhi’  ray  parameter  />  =  tit  tlA,  which 
is  a  constant  value  for  any  individual  ray  along 
its  entire  path.  The  ray  parameter  /»  or  tit  tlA 
i>  the  inverse  apparent  velocity,  or  the  slowness 


(22) 


(23) 


«A.|  2  ~  («2  +  Ml).  2 


(20) 
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nf  the  trace  nf  tlic  ray  projected  to  the  earths 
surface  at  r  =  R  and  is  defined  as 

dt/dA  =  111.1 95u(r,  6,  X)  sin  ir/R  (24) 

with  i It/dA  measured  in  seconds  per  degree,  u  in 
seconds  per  kilometer,  and  r  and  R  in  kilometers. 
It  is  a  characteristic  feature  that  this  ‘ray 
parameter'  is  no  longer  a  constant  value  along 
each  ray  if  the  velocity  is  a  function  of  all  3 
coordinates,  r  =  v  (r,  0,  A).  Thus  monitoring  n( 
dt>dA  along  the  ray  during  the  computation 
can  he  used  to  indicate  whether  the  ray  propa¬ 
gates  through  a  lateral  heterogeneity  or  not, 
depending  on  whether  dt/dA  varies  or  remains 
constant.  The  consequences  of  a  variable  dt  dA 
for  any  ray  passing  through  a  lateral  hetero¬ 
geneity  are  extremely  important  for  any  study 
using  dt  I dA  data  measured  at  large  seismic 
arrays.  The  common  hut  false  assumption  that 
r  varies  only  with  depth  can  yield  unrealistic 
or  at  least  inaceurate  results  for  the  velocity 
structure.  This  important  subject  is  further 
investigated  in  a  separate  paper. 

Before  applying  the  rav-tracing  method  to 
interpret  observed  patterns  of  travel-time  resi¬ 
duals,  an  important  characteristic  of  computed 
residuals  must  be  recognized.  Computed  resid¬ 
uals  det>end  essentially  on  the  geometry  and 
magnitude  of  lateral  heterogeneities  and  very 
little  or  not  at  all  on  the  specific  enrth  model 
used  for  reference.  Hence  it  is  irrelevant  to  the 
computed  residuals  if  a  velocity  distribution  of 
Jeffreys  and  Bullen,  Herrin  and  his  associates, 
or  any  other  reasonable  earth  model  is  used  as 
standard  for  the  laterally  homogeneous  j»art> 
of  the  earth's  mantle.  The  standard  model 
chosen  is  not  essential  Itecauso  the  travel  times 
for  i  he  standard  earth  with  superimposed 
heterogeneities  are  compared  to  travel  times  of 
the  standard  earth  without  heterogeneities  Thus 
the  travel  times  for  the  standard  earth  cancel 
when  the  residuals  are  computed 

The  situation  is  different  for  observed  resi¬ 
dual'  <  Hiserved  residual'  de|«ettd  on  the  standard 
model  i ims I  Int.-iu'c  there  is  usually  Mime  un¬ 
known  differenee  lielween  the  assumed  standard 
earth  nid  the  real  average  earth  representjng 
the  teeloniealK  Itoliaclive  parts  of  the  mantle 
Thu*,  when  oliM-ned  re-idnal*  are  obtained,  the 
travel  i itue'  for  the  actual  (bin  normals  |virts 
of  th«  earth  do  not  cancel  with  the  travel  lime* 


for  the  adopted  standard  model.  Hence  the 
observed  residuals  are  not  purely  generated  by 
lateral  heterogeneities  but  are  somewhat  falsi- 
fien  by  a  systematic  difference  between  the 
actual  and  the  adopted  standard  earth. 

Only  the  so-called  ‘relative  residuals’  as  em¬ 
ployed  by  Mitronovax  [1060]  in  his  study  of 
the  Fiji-Tonga  area  do  not  depend  on  any 
assumed  velocity  standard.  We  shall  refer  to 
these  data  in  the  following  section. 

Application  to  Tonoa-Kermadec  Arc 

To  illustrate  the  application  of  the  ray¬ 
tracing  method  to  seismic  travel-time  studies, 
an  example  is  presented  for  a  geographic  region 
where  the  existence  of  a  lateral  heterogeneity  is 
well  established.  The  Tonga-Kermadec  arc  is 
one  of  the  most  thoroughly  studied  of  the 
active  island  arcs.  The  spatial  distribution  of 
the  seism’-ity  is  described  by  Sykex  [1966]  and 
Sykex  ct  al.  [I960];  the  seismic  focal  mecha¬ 
nisms  and  their  tectonic  implications  for  that 
region  are  discussed  by  hark x  et  al.  [1969]; 
wave  propagation,  absorption  phenomena,  and 
travel-time  anomalies  were  investigated  bv 
Oliver  and  Ixackx  [1907],  Mitronovax  [1969], 
and  Mitronovax  ct  al.  [1969], 

To  illustrate  ihc  effect  of  lateral  seismic 
heterogeneities  of  the  Tonga-Kermadec  arc  on 
travel  times,  computed  arrival  times  of  traced 
rays  are  compared  to  P  travel  times  observed 
at  local  stations  from  deep  sources.  P  arrivals 
at  teleseismic  distances  from  shallow  events  in 
the  Fiji-Tonga  region  are  also  considered.  All 
live  ol)sorva‘ional  data  used  for  comparison  are 
taken  from  the  study  of  Mitronovax  [1969], 
The  velocity  model  adopted  for  the  calculation 
is  shown  in  the  lower  part  of  Figure  3.  It  illus- 
i.ates  a  .i-rtical  section  through  the  island  arc 
with  the  plate  dipping  into  the  mantle  The 
Motion  cuts  through  an  assumed  hypocenter 
located  inside  the  plate  at  4  =  20*8  and 
A  =  !79*\V  at  a  dcpvl.  of  600  km  The  profile 
strikes  N70*W  perpendicular  to  the  Tonga 
trench  The  shortest  distance  to  the  trench  is 
alioul  4V  lo  ihe  K8K  from  the  assumed  epi¬ 
center  The  geometrv  of  the  descending  slab 
wa*  simplified  compared  to  that  outlined  in 
Figure  I  lo  make  II  compatible  with  the  slab 
geomcirv  u«ed  bv  Mitrnnnrax  1 1**69]  t’ndnla- 
H  ut'  of  I  lie  sri'mic  lone  l«Jow  |00  km  were 
removed,  yielding  a  planar  dab  dipping  at 
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Fir.  3.  (Httlltm)  Vertical  section  through  the  Tonga  art*  (*crj>emhenlnr  to  the  strike  of  tli«* 
an*  The  geometry  of  the  lithospheric  plate  lielow  h  r  100  kin  i*  simplified  to  a  planar  slab 
Xitmlicr*  on  the  traced  rav*  (solid  lines)  arc  take-off  angles  i  for  a  source  at  a  depth  at  000 
km.  (Top)  P  travel  limes  versus  epieentral  distance  for  vaiion*  cases  Tlie  solid  line  ami  o|s*n 
circles  indicate  a  normal  mantle  with  a  velocity  distribution  according  to  //cm*  rt  n I  IIMMI 
without  a  dipping  lithosphere,  the  dotted  line  and  inangles  show  a  slab  wnh  a  velocity  }'■ 
higher  than  that  of  the  regular  mantle;  the  dashed  line  and  squares  indicate  a  J',  higher  /' 
velocity.  N’uniliers  '.long  the  curves  indicate  the  take-off  angle  i  for  several  rays  at  the  Maine 
Solid  circles  won  error  bars  ami  sialion  lalieb  are  olwerved  travel  limes  from  o* 

IIMRl. 

ilsatt  S.V  for  tlepth*  of  IflO  to  7(11  km.  In  tlie  wa*  a«*nnicd  in*tdc  the  dab  m  om*  ca«e  and  a 

n|'|«r  Hll  km  tlie  dip  of  tlie  plate  ilerre.ases  V;  higher  velocity  in  another  c,a«e  In  this 

to  iImhii  .Ht'  near  the  trench  The  horizontal  tn««lel  mi  low-vrlnrtty  zone  was  axunwd  ls>- 

dimen*ion  of  the  plate  was  chosen  to  lie  l<>)  neath  ami  adjacent  to  the  volcanic  ridge  Tlie 

km  The  f*-\eloritv  di*trilnit ion  of  //.  rr<«  up|s-r  part  of  Figure  .1  «howrs  the  rofrr«pnnding 
•  I  ol  |  pw»v|  »-.t«  adopted  for  tlie  mantle  mil-  ;rave|-iiim*  curve*  calculated  with  the  rav  • 
-ide  •  he  slab,  wink*  a  7^,  higher  P  vclociiv  treeing  program  The  miIkI  line  through  oj»*n 
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circles  i>  the  regular  Herrin  travel  time  with 
no  dipping  plate  present.  The  dotted  line 
through  open  triangles  corresponds  to  P  arrivals 
with  a  H**,  higher  velocity  inside  the  plate  and 
the  dashed  line  through  open  squares  to  P 
arrivals  with  a  ~r',  higher  velocity.  The  com¬ 
puted  travel  times  show  clearly  that  the  P 
waves  arrive  increasingly  earlier  (compared  to 
the  Herrin  travel  times)  the  further  the  rays 
propagate  through  the  fast,  dipping  .ithospheric 
plate  Thus  rays  emerging  at  larg  >r  distances 
close  to  the  trench  are  associated  with  larger 
negative  residuals  (up  to  —  ft  sec)  than  those 
farther  away  from  the  trench  and  closer  to  the 
epicenter. 

Superimposed  on  the  calculated  travel  time* 
are  ooserved  travel  times  from  5  stations  in 
the  Tonga  island  are.  To  obtain  tlie  'observed' 
travel  times  for  those  stations,  ‘relative  resid¬ 
uals'  determined  by  MUrttnova*  [I960]  were 
subtracted  from  the  Herrin  travel  times  at 
corre.-|>ondiitg  epicentral  distances  Mitmnovas 
attempted,  first,  to  eliminate  the  source  bias  for 
all  events  studied,  u«iug  only  those  near  and 
distant  stations  for  the  relocation  for  which  the 
seismic  rays  bypass  major  lietcrogeneities  in 
the  upper  7H0  km  of  the  mantle  in  the  Fiji- 
Tonga  region  Second.  Ik*  obtained  'relative 
re*idual*‘  for  local  station*  in  the  Fijt-Tonga 
region  The  relative  residual*  result  from  a 
comparison  of  arrival  time*  at  similar  epicentral 
distance*  itt  the  Tonga  and  in  the  Fiji  Islands 
for  ray*  that  pa**  through  or  l»y|*a**,  respec- 
tivolv.  the  dipping  plate  The*e  relative  residuals 
for  the  local  *i  tinn*  are  to  a  fir*t  approxima¬ 
tion  independent  of  the  standard  travel  time* 
n«cd  |cg.  Jrfjn  ij*  ntuf  Hullrn,  1040.  |!ttx. 
Ih  rrui  it  nl .  |1tt'»>|  and  lienee  are,  to  a  first 
order,  ill* leoeodctll  of  tin1  exact  knowledge  of 
the  e\ad  origin  tunc  of  tlie  earthquake  The 
residnab  n-*d  in  Figure  .1  to  olitain  < lie  *ol  served 
lr.a\el  lime’  are  die  average  vahle*  obtained  b\ 

. . .  i-  from  all  available  P  data  ai  the  5 

local  *1  1 1  loll-  f«  *r  ■<*  carefully  relocated  eartll- 
•|ilikc*  Willi  hi  .H cr.ige  focal  depth  of  a  la  wit 
Nil  kill 

Tin  eniiip'«risoii  iif  observed  ami  caloilateil  P 

•  roil  nine-  show  lli.il  I  lie  xc|ocil\  ltl*lde  the 

•  lab  i-  nil  i he  average  al*iin  •>',  higher  than 
l|w  III. (Ill h  vclorilie*  mil-lde  the  plate  Note 
ilid  iln«  r< -oh  i*.  io  a  tir-i  or>|er.  ind*  |»  ndeiii 
ot  ib.  m-*  ■■(  a  s|wcific  standard  nusk  I  tin  thi» 


case  Herrin  et  nl.) .  There  is  some  indicati  n 
that  this  velocity  contrast  increases  from  aliout 
«*irJ  to  alxmt  7^  as  t ho  emerging  rays  approach 
the  irench.  If  this  trend  is  confirmed  in  future 
studies,  three  different  explanations  should  lie 
considered.  First,  the  change  in  velocity  con¬ 
trast  is  real  and  such  that  it  decreases  with 
depth  inside  the  plate  liecau.se  of  a  gradual 
heating  of  the  cooler  plate  as  it  moves  down¬ 
ward  into  tlte  warmer  mantle;  hence  rays  pene¬ 
trating  through  only  the  lower  part  of  the  plate 
are  slower  compared  with  those  penetrating 
through  the  upper  part  with  the  higher  velocity 
contrast.  Second,  rays  emerging  at  smaller  epi- 
centr.il  distances  inav  propagate  a  longer  dis¬ 
tance  through  a  Imdy  of  lower  vel.icittcs  located 
near  and  lieneath  the  volcanic  licit,  as  indicated 
in  Figure  2  Or  third,  the  geometry  i«  more 
complex.  A  low-velocity  region  was  rat  taken 
into  acconii.’  for  the  model  calculations  shown 
in  Figure  .1.  If  tlie  low-velocity  rerion  exists, 
it  will  increase  the  estimate  for  the  average 
velocity  contrast  of  the  slab  to  a  higher  value 
of  almitt  7^ . 

Tlie  method  of  relative  residuals  is  es|iecially 
suitable  for  deep  events  if  local  stations  are 
available.  A*  shown  by  Mitrnnova*  [1960],  tlie 
method  of  relative  residuals  i*  not  applicable 
for  shallow  source*  (occurring  in  the  upper  j»art 
.if  the  plate)  with  Nation*  at  teleseismic  dis¬ 
tance*.  Titus  residuals  at  t«  !c«ri*mic  distances 
for  shallow  earthquake*  occurring  in  the  upper 
|ia»1  of  tlie  plate  are  dependent  on  tlie  earth 
model  it«ed  for  reference  Figure  -I  illustrates 
a  tvptcal  ca«e  showing  the  geometry  of  the 
plate,  tlie  source  location  relative  to  tlie  plate, 
and  two  selected  ray*  traced  through  the  model 
Many  rav*  striking  only  |ieqiendicular  to  the 
are  were  traced  to  teh-*ei«mic  distance*  and 
tlieir  P  residual*  wore  calculated  by  n*ing  a 
Herrin  model  for  tin-  mantle  outside  tlie  slab 
and  a  7'I  higher  velocity  in«idc  the  *lah  The 
calculate*!  r*-«t(l»nl*  obtame*l  from  rav  facing 
are  plot ual  in  Figure  .I  (*oImI  triangi.  I  a-  a 
f unci  Hill  of  epH-entral  distance  Note  again  that 
llie«c  miM/iKfo/  rc*nlual«  are  not  affecti-d  b\ 
i Ih-  n*e  of  Herrin's  model  For  comparison,  the 
*doerve*|  resHlnal*  of  ail  earthquake  ih>if  i)h 
Tonga  Island*  rcjiortid  by  I/ifronoi  tv  [  | •  | 
are  plot  let  I  t«oli(|  line  through  o|ien  circle- 1 
TI»C  loealH.n  of  till*  enrlh(|U.'lke  W||||  res|*-ct  to 
•  lie  -lab  i-  verv  similar  io  tint  *hown  in  >hc 
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model  of  Figure  4  Tints  the  two  wts  of  rcstd- 
iiaU,  the  nlniLilrd  and  the  nltscrvcd  on**,  ran 
hr  pm|icrly  iwd  for  mmpari««in.  Tlir  eliarar- 
tcrisiir  parameters  of  thi**  earthquake  arr: 
AugiM  12.  1067,  Oiili  39m  44.1s,  h  =  134  km. 
4>  -  24.Hfl‘S.  A  =  177.01 ‘W,  m.  =  .Vs  Tin* 
location  and  origin  lime  were  derived  by 
Mitnmovas  hv  using  data  from  tlir  local  *ta- 
t  tons  onl>  Tlir  residual*  wrrr  obtained  oiilv 
for  rays  traversing  through  tho  plate  in  a  \VX\V 
directum  to  lelcseismie  distance*.  Jeffreys- Bullrn 
( .I-H 1  travel  times  were  used  for  reference. 
These  J-B  residuals  were  averaged  in  10*  inter¬ 
vals  of  epieentral  dist.me* .  The  oltservcd  J-B 
ami  the  calculated  residuals  shown  in  Fig¬ 
ure  A  have  in  common  that  they  incrcat**  with 
epieentral  distance  as  the  rays  p<ii|tngate 
along  an  increasingly  longer  path  inside  the 
slab  llavs  that  travel  the  longest  |m**ililc  dis¬ 
tance  insole  the  slab  emerge  at  epieentral  dis¬ 
tances  approximately  between  AO*  and  60*  and 
show,  as  expected,  the  largest  negative  resid¬ 
uals  Beyond  HO*,  the  residuals  ilccrrase  as  the 
rays  leave  the  slab  on  its  lower  boundary.  It  is 
evi«lent  from  Figure  A  that  the  calculated  resid¬ 


ual-  do  no)  decrease  rapidly  enough  in  that 
di-t.iiuc  range  Tin-  l»  ls*can*e  I  1m*  awmrned 
horizontal  dmicu-ioti  (JtHI  knit  for  that  nnslel 
i-  hm  large  Ibalncing  tin-  thlckne—  of  the  plate 
Ito  alsnii  1 4 M I  km  a-  in  the  pretmn-  rmsleli 
would  n**uh  in  a  rum*  for  tlw  calculated  re»id- 
•lal-  that  ha-  approMliialelx  tile  -aim*  *lia|s*  a 
the  curie  for  Hie  ol»>cr\cd  re»ld'lal»  I’afl  of 
llw  large  tiegal n e  residuals  at  epirenler  ill- 
lane**-  iwar  141*  could  Is*  interpreted  alien  l- 
tiieh  bv  tlw  retain  el\  large  .1-14  irai'4  linw- 
(approximately  2  -**c  longer  than,  eg  .  Herrin 
travel  i ime- 1  In  onr  opinion  tin-  i-,  however, 
more  likely  an  indication  that  the  travel  time- 
of  Herrin  et  al  an*  somewhat  t«si  short  in  ilu- 
*li»iaiice  ring**,  since  ilwy  do  not  yield  tlie  largr 
negative  n --idu.il-  i*x|iected  from  the  stu*l\  of 
ileep  events  in  that  ngmti  It  i>  not  our  inten¬ 
tion  to  determine  in  this  |ia|ier  tlw  *ha|ie  ami 
siz**  of  tlie  flipping  slab  in  tlw  Tonga  region 
bv  I  mug  the  nnslel  «o  that  oImtiisI  anil  cal¬ 
culated  re-nlnal-  match  m  an  optimum  sense 
However,  it  is  clear  that  by  systematic  varia¬ 
tion  of  tlw  six  nnslel  iiarameter*  and.  if  twee*- 
sarv.  of  ilw  sha|«*  of  tlw  *ei*mw  zmw  one  ran 


fig  4  A  nnslel  of  the  slab  with  a  -ofiro  al  dipih  k  100  km  Tin-  lists  h  I  i-  n->  ,|  (nr 

inning  ravs  III.  an<l  •  aknl.i'ing  /’  rr«ldi|al«  at.  If  k  «•  1*011*  •l»*».an«fs  ISfi.  Ilf  i«n  to  ft  in  « 

witii  islcr.al  r>(  rail  SKIS  .an  s|io«n  Th<-  «aloilit«*l  /'  r>-l>l<lfl-  lor  Mill  iimhIs  I  ir<  (.kilt.. I  m 
Figure  5 


SKISMIC  RAY  TRACING 


Mi87 


Fig  5  Cnin|<in<nn  of  olumnl  /*  mndiial*  I.Wtlmwii'iw.  IMOI  nod  iheorrinal  re»ninaU  r»|* 
i  illaled  |»y  ray  I ruing  for  i Ik  mmirl  >hi>«n  in  F’lgUfe  I  F'uf  ilelniU  m-r  lr\l 


find  a  nmilnn.il Hin  of  |«arametrr»  lh.ii  will  yield 
.1  minimum  for  the  »i.ind.ird  deviation  Iwtwern 
oliwerved  .iihI  ralriiLited  re*idital»  Thi*  method 
of  trial  .iml  rrmr  may.  Ihii  ilom  not  nermaanh . 
converge  toward  the  arMiui  vloriiy  atntetitrr, 
I  n  Hu  IK  ImMii***  of  the  re»tnrtmn«  inherent  in 
i he  tnmlrl  Literal  hcicrogmcttie*  nmr  the 
nrrmr  or  .imwlterr  along  tin-  mi  of  the  ray 
|*alh.  and  |mM|li|r  difference*  Irtnwn  the 
ai|o|*iii|  *i.ntdard  niodi’l  and  tin*  ariu.il  earth. 
mm|*liraie  iIh-  proMem  ami  introduce  non* 
nni<i<K  iH  *•  Further  uneviKatmiH  ate  ncrc**arv 
to  nm«nlrr  more  rtgorou*  inversion  technique*, 
eg.  iIhw  i|r\  rln|Mi|  l<\  Hack  u»  am/  ( iilbrrl 

|  1W7.  I'«»«»| 

Tin  rr  I*  a  |H  -lull  of  almill  I  A  *cr  I*- 
IWiHl  lie  Iilw|\n|  llld  I  In  r.ilrillaled  r»'“l«l- 
oal*  m  Figure  A  Tin*  nff-ci  I-  |«nm.inK  due 
in  *  In-  unnn.iitiM  in  iIh’  origin  tutu-  of  i  !*•' 
i  ,iftln|i|  <k<  iihI  ikI*  oil  i In-  *|irr|fir  t|*«-  of 
i  In  I  I  *  •  •!•!*-  i*  i  In  *  •  r  <\ •■l-iiiit*'  -ti  initial  for 
r<  l'«  linn  oi  iIh  in  in  iihI  ilrirmiin.il  n«n  of 
*•  iiioii  n  *k|ii  iI*  Mon  n -li.il nifornriiioti  em 
In  i  \i rii  iiil  i r < hi i  n  -kI-hI*  of  large  nuele.ir 

•  \  |  d- ’*  i*  in  *  lor  w  1 1  h  1 1  i  Ik  one  iii  nnn  md  Im-.i- 

*  ion  ir>  known  1 1  n  ir  iii-K  TIk’  /’  r*'*idu.il« 


of  the  nnrlear  underground  c\|i|n*inn  Lmgdnit 
will  lir  iimiI  in  a  later  |ia|ier  to  invent «ate 
lateral  hrirrogenrUir*  awortatr*'.  wrilh  the  Alett* 
i kin  i*l.i ml  arr 

Sovr  Sfct«\HH^«.U  AL  lMn.ir*TtoNa 

A  lw«ie  numertr.il  *rhrmr  ir  |irr*ented  in  tin* 
|ia|irr  for  tearing  *et*mie  tav*  through  a  Liter* 
alK  lietengpiteiHi*  earth  Willi  ihn  new-  method 
it  |.  |H«**ili|r  for  i lie  fir«t  time  in  r.ilriil.ite  on  a 
worldwnic  *eali'  areitratr  thenmieal  muditaU  of 
in  Iter  /*  or  S  travel  lime*  r.m*ei|  K\  a  glolial 
|>alli'tli  of  lateral  lielrfrigeneit le*  in  tile  rniM 
■Hid  n|i|ier  manile  Tile  gnwnctnr  •liape  of  I  lie 
lerfonic  mill*  and  lielemgetteiHI*  Imdir*  ran  lie 
derm-d  for  a  hr*i  trial  model  from  \ anou* 
Miurre*  of  mfomt.il  Hill.  eg.  from  «rl«mie|tV.  the 
gi  gnplm  di*inlHiltoit  of  i|rr|t**r.i  innelic*. 
ari|\e  \ oli  HKK'*,  mnl*orcanir  mlgf*.  and  otlier 
gro*.  lerlolilT  altil  mor|ihol-igieal  feililO-*  Hv 
*\ »leni.llie  ro|il|<:in*o(l  of  I  hmTct  tea  I  re*ldoil* 
iihI  oliM-rxial  r>-*idii.il*  from  nuclear  iA|i|o*nin* 
1*  Will  I*  earl Ih|U.iLi*  in  \anoil*  lerlonic.illv 
a  ini'  region*  around  die  world.  iIh  |iar.inieier* 
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lirrent  seismological  studies  suggest  Intend  inhomogonoitios  in  /'  and  ,S  velocities  of  the 
1,1  dull,  are  associated  with  slabs  of  mobile  lithosphere  descending  into  the  mantle  lie- 
lie, ilh  island  ares.  In  special  cases,  travel  times  of  1‘  traversing  sueli  zones  ran  differ  l>v  as 
much  as  5  sec  ami  of  .S'  by  up  to  10  sec  from  standard  travel  times.  In  addition,  such  zones 
are  characterized  by  relatively  low  attenuation  of  .S'-wavc  energy  compared  with  high  atlenua- 
Hon  in  a  broad  zone  oil  the  landward  side  of  the  active  volcanoes.  To  explain  the  observed 
anomalous  travel  times  and  attenuation  phenomena,  it  is  necessary  to  trace  the  path  of  body 
waves  through  laterally  heterogeneous  earth  models.  The  techniipie  of  rav  tracing  developed 
here  uses  Fermat's  principle  to  obtain  the  differential  equation  of  a  ray  in  spherical  coordi¬ 
nates.  I  lie  position,  direction,  and  travel  time  of  the  seismie  'wave  front  at  any  point  along 
the  curved  ray  path  are  obtained  hv  mmierical  integration  of  the  differential  equation  for  an 
assumed  three-dimensional,  continuous  velocity  distribution.  The  proh'eni  of  representing  a 
realistic  three-dimensional  velocity  structure  in  the  earth  is  solved  in  a  way  that  is  especially 
suitable  for  use  on  computers.  Some  examples  for  rays  traversing  an  island-arc  structure  are 
presented.  The  implications  of  this  method  of  tracing  rays  in  a  laterally  heterogeneous  earth 
are  discus-,  , |  with  respect  to  seismic  travel-time  studies,  interpretation  of  residuals  in  terms 
ul  b  i  iot.:  'i  rogeneities,  source  bias,  and  the  precise  location  of  earthquakes  mid  nuclear 

i  xplosioi  /a  Measurements  from  large  seismic  arrays  and  their  inversion  to  obtain  de¬ 
tails  ol  i!  a'iiy  strm  rime  in  tlx*  up|>er  mantle  are  also  discussed. 
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